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We study the dynamics of a Bose-Einstein condensate (BEC) in a one dimensional optical lattice 
in the limit of weak atom-atom interactions by incorporating quantum fluctuations. The pulsat- 
ing dynamical instability manifests itself in the time evolution in which atoms periodically collect 
themselves into a pulse and subsequently disperse back into the initial homogeneous state. We 
take into account the quantum fluctuations within truncated Wigner approximation and observe 
that the quasiperiodic behavior still persists for single realizations which may represent the typical 
experimental outcome. The quantum mechanical ensemble averages of the wave functions shows a 
damping in the pulsating event. The fluctuations become more prominent for smaller atom numbers. 
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I. INTRODUCTION 

The superfluidity of a Bose-Einstein condensate (BEC) 
in an optical lattice has been drawing a considerable at- 
tention in last several years [l[ . As is well known, super- 
flow of the BEC in free space suffers from an instability 
when the center of mass velocity reaches a critical value. 
Such an instability, known as Landau or energetic insta- 
bility, exist when the supcrfluid flow is not at a local 
minimum of energy and the system lowers its energy by 
emitting phonons [2|. In an optical lattice, in addition to 
the energetic instability, the BEC may also exhibit dy- 
namical or modulational instabilities which have been a 
subject of active experimental and theoretical research in 
recent years]! 1, i, 1 0, I, ©, EE El El El, EI EE El, 
E3j EM El) l20j - When the system is in the dynamically 
unstable regime, small perturbations grow exponentially 
in time resulting in an irregular dynamics, loss of coher- 
ence or an abrupt stop of the transport of the atom cloud 

In this paper we study dynamical instabilities of atoms 
in an optical lattice for the case of weak atom-atom in- 
teractions and also taking into account quantum fluc- 
tuations of atoms. We recently reported [2(| that, by 
appropriately selecting the strength of the interactions, 
the corresponding classical system may exhibit a pul- 
sating dynamical instability in which the atoms nearly 
periodically collect to a peak in lattice occupation num- 
bers, and subsequently disperse back to (very close to) 
the initial unstable state. This is different from the con- 
ventional view, valid at strong interatomic interactions, 
that dynamical instabilities for BECs in optical lattices 
are associated with irregular dynamics. When we incor- 
porate quantum fluctuations of atoms using stochastic 
phase-space methods, the quasiperiodic behavior is still 
observable in individual stochastic realizations that rep- 
resent typical individual experimental realizations. As 
the pulsating solitons in each realization appear at dif- 
ferent lattice sites due to quantum effects, the quantum 



mechanical ensemble averages of the wavefunction revival 
become progressively weaker when the effective interac- 
tion strength is increased. Other ensemble averages, such 
as the pulsation amplitude, can still provide information 
about the quantum soliton. 

We consider a stationary superfluid flow of a BEC in 
an optical lattice with a large enough flow momentum 
that triggers the dynamical instability of the correspond- 
ing classical nonlinear system. In a quantum system the 
corresponding sharp transition to the dynamically un- 
stable regime is smeared out, typically resulting in a pro- 
gressively increasing dissipation in th e dy namics close to 
the classical onset of the instability [15| ■ We provide a 
qualitative explanation of the pulsating phenomenon by 
studying the dynamics of an integrable double-well sys- 
tem. Although the instability is a result of the interplay 
between the lattice discreteness and the nonlinearity that 
makes the lattice non-integrable, the dynamics of the lat- 
tice with many sites is approximately as if the system is 
integrable. Related classical pulsations starting from al- 
ready compressed atom distribution in a lattice have been 
discussed in [l8| within the frame-work of the nonpoly- 
nomial Schrodingcr equation. 

The pulsating instability manifests in the dynamical 
regime where the nonlinearity is weak. In the mean-field 
description, the size of the nonlinearity is proportional 
to the atom-atom interactions and the total number of 
atoms present in the system. As the number of atoms 
gets small the mean-field description may breakdown as 
the relative fluctuations in the system amplifies, and the 
quantum treatment is inevitable. We would like to know 
how the quantum effects smear out the pulsating mech- 
anism as we reduce the number of atoms in the system. 
At the simplest level, we study the quantum dynamics 
of the pulsating instability using the quantum distribu- 
tion function, in particular, the Wigner function method. 
The Wigner method simulates the quantum mechanical 
system in classical stochastic process where the quantum 
fluctuation is included in the initial state. In the case 
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of BEC it gives the time evolution of the whole mat- 
ter held including both condensate and non-condensate 
atoms, and allows the scattering between them, which is 
absent in the classical GP description. 

In Sec. II we formulate the theoretical model, mainly 
the Gross-Pitaevskii equation (GPE) [23] and its dis- 
crete variant, the discrete nonlinear Schrodinger equa- 
tion (DNLSE) dHI. We use linear stability analysis to 
hnd the region of interaction strengths and flow quasi- 
momenta where the system develops instability. As in 
nonlinear dynamics, following [30j . we verify the exis- 
tence of the localized soliton solution in the forbidden 
gap of the linear spectrum. In Sec. Ill we investigate the 
time evolution of the DNLSE within classical mean-field 
theory. Although the system initially develops instability 
the time evolution shows a regular dynamics whereupon 
the atoms periodically collect themselves into a pulse and 
disperse back into the unstable state. 

In Sec. IV we review the well understood double well 
system and argue that the dynamical behavior of the 
multi-site system is analogous to the two-site system, at 
least in the limit of weak nonlinearity. In Sec. V we study 
the dynamics beyond the classical mean field theory us- 
ing truncated Wigner approximation (TWA); a phase- 
space method that approximately solves the dynamics of 
a quantum system by means of stochastic initial config- 
uration. We then compare various physical properties 
such as the number fluctuations and the overlaps of the 
state of the system in single realizations with an ensemble 
averages. Quantum dynamics significantly modifies the 
classical picture as the number of particles gets small. 
We observe the damping in the pulsating phenomenon 
when we average over many stochastic trajectories. 

II. THEORETICAL MODEL: DNLSE 

At absolute zero temperature the dynamics of the BEC 
atoms in an optical lattice can be modeled by the mean 
field Gross-Pitaevskii equation 0, HtJ 

where ^(x, t) is a wave function corresponding to the 
bosonic field operator such that ^(x)! 2 = n(x), the atom 
density. The coupling constant g^D is related to the 
scattering length through g 3 o = 47r ^ n a " , where a s and 
m are the s-wave scattering length and atomic mass re- 
spectively. The positive and negative scattering lengths 
respectively correspond to the repulsive and attractive 
atom-atom interactions. The Eq. ([T]) is an approximate 
description of an assembly of a large number of bosonic 
atoms that are in the same quantum mechanical state. 
We consider the external potential, V(x) of the form 




Here Vq and cZl(= A/2) are respectively the depth and 
the periodicity of the optical lattice. If the harmonic 
confinement is much stronger in the transverse than in 
the longitudinal direction (ui± ^> ui x ) the GPE can be 
transformed into a one-dimensional form 



with an effective atom-atom interactions g — 2a s hu)±. 
When the depth of the optical lattice is much larger than 
the chemical potential of the atoms, one can employ the 
tight-binding approximation. By expressing the conden- 
sate wave function ip(x) as a superposition of the Wan- 
nier functions localized within each potential well of the 
lattice, one can obtain the tight-binding version of the 
GPE known as the discrete nonlinear Schrodinger equa- 
tion (DNLSE) @: 

d 2 

l "a7^« = ~ J Wn+l + Vpn-l) + (Vn + X Wn\ Wn- ( 4 ) 



The parameters J and V n respectively characterize 
the tunneling rate and the external trapping potential, 
whereas x is proportional to the atom-atom interactions. 
It is convenient to scale the wave function and interac- 
tion parameter as ip n = \^rtp n and x = NtX> with Nt 
being the total atom number, so that the wave function is 
properly normalized to one. Unless it is explicitly stated 
otherwise, we assume here that J > and the atom-atom 
interaction is repulsive, x > 0. 

In the absence of external potential and nonlinearity 
the Eq. ^ may be solved with a plane-wave ansatz 
V-Vi(£) = -i=e J ( pn ~"( p '*' , giving the dispersion relation 
u){p) = — 2Jcosp. The periodic boundary conditions 
quantize the quasimomenta p — 2ttP/N, N being the 
number of lattice sites, and P is an integer that may be 
chosen to lie in the interval [-y, y)- For notational con- 
venience we always take the number of lattice sites N to 
be even. 

When the interaction is switched on, the constant- 
amplitude plane waves are still solutions to Eq. Q but 
with a modified dispersion relation u(p) — — 2 Jcosp+-£. 
Besides these extended- wave solutions, the nonlinear sys- 
tem also admits solutions that are localized in space [311 ] . 
These solutions, so called the gap solitons, usually have 
an energy that lies outside of the linear band spectrum. 
In the continuum model solitonic solutions of the nonlin- 
ear Shrbdinger equation can be found in closed form by 
using the inverse scattering method [32]. However, the 
discrete system has fewer constants of the motion and is 
not integrable as such, so that one has to rely on numer- 
ical techniques. 
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A. Hamiltonian and its symplectic nature 

The Hamiltonian corresponding to the equation of mo- 
tion Eq. (HJ) is given by 



H = Y, \ -J(M; i+1 +h.c.) + V n \yj n \ 2 + ^„ 



(5) 



where ip n and itp^ are canonically conjugate variables 
that satisfy Hamilton's equations of motion 



1pn = - 



dH 



dH 



(0) 



Although we are dealing with a system with a quantum 
origin, the macroscopic wave function, nonetheless, obeys 
classical equations of motion. Since the time is cyclic in 
the Hamiltonian, the total energy is a constant of the 
motion. The normalization which is proportional to the 
total number of particles, is also a constant of the motion. 

To study the basic features of the solutions governed 
by the Hamiltonian (Eq. |(SJ)) near the edge of the lin- 
ear band spectrum, we consider Eq. Q as a map where 
the lattice indices play the role of the discrete time [30l ] . 
Without the loss of generality and for simplicity in our 
continuing discussion in the following we neglect the ef- 
fect of external trapping potential so that the system is 
translational invariant along the lattice direction. First 
we write down the stationary state solution of the Hamil- 
tonian in the form tp n (t) = ip n exp(— iut) to obtain the 
time independent equation 



-J(lp n+l +1p n -x)+X\lpn\ 



(7) 



and then separate the real and imaginary parts, Vn = 
x n + iyn, resulting an area preserving 4-dimensional real 
map M : 



j n 

Vn 



•En i 
Vn- 



vl)Xr. 
vDVn 



■ Ur. 
V„ 



(8) 



Here, for convenience, we take J = 1. Given the initial 
conditions [xq, ya,UQ,Va), one can propagate the solution 
for a given energy oj to obtain an orbit of the discrete 
lattice system. The Jacobian matrix of the map M is 
given by 



fab 
b c 
1 

Vo i 



-1 

3 -1 

3 

3 



where 



a 
b 

c 
7 



7 + 2xz„, 

nVn > 

i + ^xvl, and 
x{%l + vl) - 



Since the determinant of the Jacobian matrix is one, the 
map is indeed area preserving [33j . The fixed point of 
the map is (0,0,0,0). In order to study the stability of 
the fixed point one has to solve for the roots of the char- 
acteristic polynomial 



(A(w + A) + 1) 2 = 0, 
which gives the corresponding eigenvalues 



A j 



-oj ± {oj 2 - 4) 



1/2 



(9) 



(10) 



Since A + A_ = 1, the roots are reciprocal of each other in- 
dicating the symplectic nature of the Hamiltonian. There 
are three possibilities: 

(a) —2 < oj < 2; all the roots are complex with magni- 
tude one; 

(b) oj > 2; the roots are real and negative; 

(c) oj < —2; the roots are real and positive. 

In case (a) the periodic orbit corresponding to the fixed 
point (0, 0, 0, 0) is elliptical and is stable. In case (b) 
and (c) the periodic orbits are hyperbolic and is unsta- 
ble. In Fig. [1] we have plotted an orbit of the map M 
for oj = 1.95 (top) and oj = 2.05 (bottom). It is clearly 
seen that the elliptical fixed point (0, 0, 0, 0) loses its sta- 
bility and turns into an unstable hyperbolic point with 
the onset of period doubling bifurcation when oj passes 
through the critical value two. This hyperbolic fixed 
point lies on a homoclinic orbit that corresponds to spa- 
tially localized soliton solution. Furthermore, it should 
be noted that Eq. J7]) is invariant under transformation 
ip n — > (— l) n ip n , x ~ > ~Xi u ~ * ~ w j every solution in 
the positive region of the linear band spectrum has one 
to one correspondence to the negative one. We, thus, ex- 
pect soliton solutions in a lattice for both repulsive and 
attractive atom-atom interactions. However, these two 
solutions differ intrinsically in the sense that the solitons 
in the repulsive case have an alternating signs between 
adjacent lattice sites whereas in the attractive case, they 
have the same sign. 



B. Modulational Instability 

The DNLSE admits stationary solutions of the form 



~2Jcosp+ j^. 



^ e i(pn-u;(p)t) with the dispersion = 

. To study the stability of a solution we in- 
troduce an infinitesimal perturbation around the steady 
state [3, [![> 

ijj n (t) = + ue liqn - nt) + v*e~ t{qn - n "% (11) 

where q and O are the momentum and the frequency of 
the small excitation relative to the initial unperturbed 
steady state solution. After inserting Eq. (TiT]) in Eq. 
(j4|) with V n = 0, and expanding to the lowest nontrivial 
order in u and v we get the following matrix equation, 



(12) 



4 



0.06 



0.04 - \ 



0.02 



Substituting Eqs. 
obtain 



a -Mu 



(16) 



2Q-Mu-M 2 2 
HI), Ha into Eqs. flg), ffTB]) , we 
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N 



(17) 
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FIG. 1: Map orbit around an elliptical fixed point for two 
energies u) = 1.95 (top) and 2.05 (bottom). The stable el- 
liptical fixed point bifurcates leading to the period doubling 
bifurcation when the energy crosses through the critical point. 
The changing of the elliptical fixed point into hyperbolic is an 
indication of the emergence of different type of solution. 



where £ is a vector [u, v] T and M is a 2 x 2 matrix with 
elements 

Mu = 77 + 4sin f sin (§+P), 
X 22 = -^-4sin| sin(|-p), 



Mu 



-Ml, 



N' 



(13) 



The eigenvalues of M give the small-excitation frequen- 
cies, 

= 2Jsinpsinq 



±y 4Jcosp(l — cos<?)[— + Jcosp(l — cosg)],(14) 

whereas the eigenvectors give the corresponding mode 
functions: 



fl-M 2 2 



2Q-Mu - Mi 



(15) 



'4Jcos(p) sin 2 (§) -T- 



N 



2T 



(18) 



with the definition 



r = ±W4Jcosp(l 



Jcosp(l — cos q)] . (19) 



By inspection it can be easily verified that the mode 
functions u and v satisfy the normalization condition 
\u\ 2 — \v\ 2 — ±1 as long as T is real. 

The eigenvalues fi corresponding to positive normal- 
ization gives the physical small-excitation frequencies 
whereas that corresponding to negative normalization 
are unphysical. Since we are dealing with the repulsive 
atom-atom interactions, x ^ 0, it is seen that all eigen- 
values are real for \p\ < ^. However, the existence of 
complex eigenvalues cannot be ruled out in the interval 
\p\ G , 7r] , depending on the values of J, x, q and p. 
When an eigenvalue is complex, i.e., T is imaginary, small 
perturbations in the steady flow grow exponentially in 
time. In this case the norm of the eigenvector \u\ 2 — \v\ 2 
vanishes identically [4j. 

For any \p\ greater than ^ and for large N , the eigen- 
frequencies will be complex if 



A > | cos(p)| 



N 



(20) 



where Q is a non-zero integer that lies in the interval 
[-j, -y) and the rescaled interaction strengths A is de- 
fined to be 



A = A 



X_ 
2J' 



(21) 



The flow with the quasimomentum p is then said to 
be dynamically unstable in the sense that a small noise 
drives the system far away from the equilibrium state. 
It should be pointed out that the critical interaction 
strength approaches zero when the number of lattice sites 
goes to infinity, implying that any flow with \p\ > -| and a 
fixed A > will turn unstable with TV — > oo. The dynam- 
ical instability can be qualitatively understood from the 
dispersion relation of the DNLSE. For \p\ > |, equiva- 
lently when the effective mass is negative, the interaction 
shifts the frequency upward in the forbidden gap of the 
linear spectrum where the plane wave solution cannot ex- 
ist, which means that the system is unstable. Moreover, 
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FIG. 2: Coefficient of exponentiai gain G for an excitation 
mode with quasimomentum q for the interaction strengths 
A = 0.32 (a) and A = 5.0 (b). Modulational instability is only 
possible for G > 0, and excitation modes may only occur for 
nonzero integer values Q of q. 

the imaginary part of the complex eigenfrequency as well 
as the corresponding eigenvectors are the same for q and 
— q (see Eqs. (fl7|) . (fl8|) ), which indicates that these two 
modes are equivalent as it comes to the instability. 

Fig. [2] shows the gain curve, G — min[0, |Im £l(p, q)\], 
for two interaction strengths (a) A = 0.32, and (b) A = 
5.0 in a lattice of 32 sites with p = tt for — oo < q < oo. 
The figure (a) reveals a single pair of sidebands with one 
unstable mode, whereas the figure (b) shows four pairs 
for four unstable modes. 



III. TIME EVOLUTION AND PULSATING 
INSTABILITY 

We carry out numerical simulations on the DNLSE to 
study the growth of the unstable mode in a lattice for 
a suitable range of interaction parameters. For a given 
number of lattice sites and the flow momentum the num- 



ber of unstable modes in the linear stability analysis de- 
pends only on the interaction parameter A. Here we focus 
only on low energy excitations in the limit of weak atom- 
atom interactions. Two numerical methods have been 
used for the time evolution, an unconditionally stable 
Crank-Nicholson type algorithm [34| and a sixth-order 
accurate FFT split operator algorithm that works in the 
same way as is discussed in [351 ) for the ordinary nonlinear 
Schrodinger equation. 

A. Single Unstable Mode 

A straightforward analysis of the eigenfrequency ex- 
pression Eq. (Tj3| suggests that the range of interaction 
strengths where the Q = 1 mode is unstable but the 
Q = 2 mode is not is given by 



|cosp| — < A<4|cosp| — . (22) 

Fig. [3] shows a typical density plot of the time evolution 
of the BEC initially prepared in the plane wave state 
at the edge of the Brillouin zone (p = tt) seeded with 
random Gaussian noise, for the number of lattice sites 
N = 32 and the interaction strength A = 0.48. This 
value of A corresponds to one unstable mode in the linear 
stability analysis. Although a tiniest amount of noise 
(either in real experiments or in numerical simulations) in 
the unstable direction triggers the instability, an external 
noise of amplitude £ = 10~ 4 is added just to speed up the 
instability. It has been tested in a number of runs that 
the time for the onset of the instability for fixed values 
of the other parameters depends logarithmically on the 
amplitude of the added noise. 

Fig. [5] depicts a snapshot of a pulse that moves during 
its formation from the initial flow state with quasimo- 
mentum p = ^j^-- In this figure we take a larger lat- 
tice with N = 128 sites and the interaction strengths is 
A = 0.25, so that there is still one and only one unstable 
mode. By virtue of the periodic boundary conditions a 
pulse that goes over the right edge will reappear at the 
left edge of the lattice. 

The pulsating behavior of the peak can also be viewed 
by plotting the fraction of the initial state ip n (0) remain- 
ing in the state of the lattice as a function of time, 

f(t) = \J2r n (o)Mt)\ 2 - (23) 

n 

In Fig. [5] we plot the overlap, f(t), as a function of time 
for the same parameters as in Fig. [3] It is revealed 
that the instability drives the system far from, and sub- 
sequently brings it back to, the original unstable steady 
state, and the process repeats. Each dip in the plot rep- 
resents formation of a pulse during the course of time. 
It is also noted that the quantity f(t) does not vanish 
all the way to zero, indicating that the pulsed state is 
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FIG. 3: Collapse and revival of the density hump in the BEC 
evolution when the lattice starts from a dynamically unstable 
state. The parameters are A = 0.48, p = tt, and N = 32. 
These parameters correspond to a single unstable mode, as 
per linear stability analysis. The instability drives the ini- 
tial homogeneous atom distributions into a density peak that 
subsequently disperses back to the initial state. 




FIG. 4: Density evolution of the BEC for the initial flow 
momentum different from tt. Here the parameters are N = 
128, A = 0.25 and p = The pulsating behavior of the 

peak still persists but the peak moves with the group velocity 
v g = sin(p). 



not orthogonal to the initial steady state. Furthermore, 
a closer inspection of this plot shows that the subsequent 
peaking events are not strictly periodic; the interval be- 
tween the dips varies slightly, implying a quasi-periodic 
phenomenon. 

By analyzing data sets of this kind a number of inter- 
esting observations emerges, (i) First, in contradiction 
to the common belief that the instability may develop 
an irregular dynamics, it causes the atoms to pile up in 
a single-peaked distribution of the occupation numbers 
|t/>„| 2 . However, upon further time evolution, the system 
returns very close to the initial unstable state, again pul- 
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FIG. 5: Fraction of the initial state f(t) in the state of the 
lattice plotted as a function of time. Here the parameters are 
N = 32, A = 0.48 and p = tt. 

sates to a peak, and so on. We have periodic peaking 
and recurrences to the unstable initial state, (ii) Second, 
the peak may occur at any lattice site. It is the random 
noise that seeds the position of the peak. In order to test 
this claim, we ran the simulation a number of times with 
everything else except the particular realization of the 
noise held unchanged, and observed that peaking occurs 
approximately at the same time but the position of the 
peak is completely random. This is in accordance with 
the theory of the translational invariance of the lattice: 
A lattice-translated pulsed solution is also a degenerate 
solution of the DNLSE and there is no preferable lattice 
site for the occurrence of the pulse, (iii) Third, the peri- 
odic recurrences and the velocity of the peak for a given 
number of lattice sites, seem to depend on the values 
of interaction strengths A and initial flow momentum p 
only. For the initial flow state p =/= tt, the pulse moves 
with the velocity that turns out to be the group velocity 
of the carrier wave, v g = sin(p). 

B. Multiple Unstable Modes 

Equation (J^H) can be generalized to obtain the condi- 
tions for Q unstable modes, 

2 2 

Q 2 |cosp| 1L < A< (Q + l) 2 |cosp| 2L. (24) 

For a lattice with a large number of sites the one-peak 
condition is highly impractical because the interaction 
strength A needs to be extremely small and the pulse re- 
vival period is long. For reasonable interaction strengths, 
Eq. suggests that there may be more than one un- 
stable mode. Fig. [6] is a typical representative of the 
dynamics of the BEC for multiple unstable modes. Here 
we take N = 128, A = 2.0, p = tt and £ = 1CT 4 . Each 
bright white spot represents a pulse. As before the right 
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FIG. 7: Time evolution of the Fourier modes tp q of the lattice. 
The parameters are N = 32, p — n and A = 0.48. The mode 
Q — corresponds to the initial steady state while modes 
Q = 1, 2, 3 . . . are the low-lying excitations. Only the modes 
for Q > are shown. 



FIG. 6: Density plot of the populations \ipn\ 2 in the case 
where there are four unstable modes. The parameters are 
N = 128, A = 2.0, p = tv and £ = 10" 4 . Lighter shading 
represents higher site populations. 



edge of the plot wraps around to the left edge by virtue 
of the periodic boundary conditions. The random noise 
seeds approximately four pulses. However, these pulses 
are not independent of each other. Presumably because 
of nonlinear mode-mode interaction, they move around, 
join and split as they collapse and revive. 



C. Evolution in Fourier space 

The recurrences observed in the peaking events in the 
DNLSE resemble the energy recurrences in the Fermi- 
Pasta-Ulam (FPU) problem [H]. The FPU model deals 
with the evolution of a lattice chain with nonlinear in- 
teractions between the nearest-neighbor atoms when ini- 
tially a single low-energy mode is excited. For a time 
scale much longer than the time period of the normal 
modes, the energy is well localized to the given excited 
mode, while the amplitudes of the higher-energy modes 
decay exponentially as a function of the energy difference 
from the initially excited mode. For a longer time scale it 
has also been noticed that recurrence of the initial excited 
mode is possible. 

The pulsating behavior of the density distribution of 
the BEC atoms in the lattice can be viewed as a similar 
recurrence phenomenon as observed in the FPU model. 
We have started with a steady state for a given flow quasi- 
momenta (p > tt/2) and a suitable nonlinear interaction 
strength to trigger the instability in the system. Ergod- 



icity immediately suggests that the energy initially fed 
into a single mode should distribute evenly between all 
Fourier modes. However, the excitation amplitudes of the 
modes other than the mode corresponding to the initial 
steady state seems to decay exponentially with the index 
Q. The energy localization to a few Fourier modes in a 
nonlinear system is not a new phenomenon (3lT | . The ex- 
istence of discrete breathers in a nonlinear lattice system 
is an example. Recently, energy localization in Fourier 
space in a so called 'q-breather' has been investigated in 

Figure [7] shows the time evolution of the Fourier modes 
ip q of DNLSE for the parameters A = 0.48, N = 32. 
Only a few components are seen to be excited, as the 
amplitudes of the higher-energy modes are suppressed 
exponentially. The dominant Fourier components are P 
and P ± 1 which implies that the excitation modes that 
go unstable will have indices Q = ±1 with respect to the 
initial steady state, as expected. 



IV. DOUBLE WELL ANALOGY 

In order to explain qualitatively the pulsating behavior 
of the density distribution of the BEC in the lattice, we 
study a coupled double-well system. Assuming ipi t 2 = 
|^i,2|e 1 ^ 1,2 , the coherent dynamics of such a system can 
be described by a pair of equations [38[ , 



z(t) = -yjl- z 2 (t) sin 0(f) 
z(t) 



(j)(t) = Az(t) + 



cos (/>(t), 



(25) 



where z = \ip2 \ 2 — |"0i| 2 an d 1 



hi are the fractional 



population imbalance and the relative phase between the 
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two wells. The normalization is |?/> 2 | 2 + \^i\ 2 = 1- The 
Hamiltonian (the total energy) in these variables be- 
comes, 



(26) 



Both the norm and the Hamiltonian are the constants 
of the motion and thus the double-well system , in prin- 
ciple, is integrable. By inspection it can be checked that 
the fixed points of Eq. (|^5|) are z — 0, <p = nir, where 
n is an integer. The potentially unstable steady state in 
the multiwell system can be translated into the two-well 
system by taking the solution z = and (f> = n. The 
behavior of the orbits near this equilibrium point can be 
examined by using linear stability analysis as before. It 
can be easily verified that the state {0, 7r} is stable for 
the values A < 1, and unstable otherwise. In Fig. [5] 
we have shown the energy contours of the two-well sys- 
tem for A = 0.5 (a) and 1.5 (b). We have drawn the <j> 
axis from to 2tt so that the potentially unstable fixed 
point (0,7r) lies at the centers of the plots. For A = 0.5 
the potentially unstable steady state is an elliptic fixed 
point and the time evolution takes the system periodi- 
cally around this point. At A = 1 the elliptic fixed point 
bifurcates, and for A = 1.5 there is a homoclinic orbit 
with the emergence of two symmetric off-centered elliptic 
fixed points. Thus, starting in the vicinity of what used 
to be the potentially unstable steady state, the system 
takes off in an unstable direction along the homoclinic 
orbit and goes around one of the bifurcated elliptic fixed 
points. 

The double- well system allows an analytic solution in a 
closed form in terms of Jacobian elliptic functions. Here 
we have expressed a solution valid in the range < z < 1, 



z(t) 



(zi - z 2 )z 4 sn(i7|fc) 2 + zi(z 2 - z 4 ) 



(zi - z 2 ) sn(t7|/c) 2 + z 2 - z 4 
where we have defined 



(27) 



zi 



Z3 = 



p + a 



p — a 



z 2 = 



z 4 = 



p — a 



p + a 



2aH - 1 



a = \/p 2 + 4q, 



1-H 2 



q = 



A 

a = — , 
2 



7 = aV ( z i - z 3 )(z 2 - z 4 ) 
and 

(zi - z 2 )(z 3 - z 4 ) 



k = 



(zi - z 3 )(z 2 - z 4 ) ' 



and H stands for the conserved value of the Hamiltonian. 
The norm which is also a conserved quantity is taken to 
be equal to one. 



1.0 



0.5 



W o.o 



-o,s 



-1 — i— i — r— i — [- 




FIG 



8: Contour plots for the two site Hamiltonian in the 
z) plane for interaction strengths A = 0.5 (a) and 1.5 (b). 




FIG. 9: Time evolution of the population imbalance z given 
by the equation Eq. (|27[) f° r the parameter A = 1.5 and the 
Hamiltonian H — 1.0. 



In Fig. [TT]we have plotted the function z(t) given by 
Eq. (|27|) for the parameters A = 1.5 and H = 1.0 such 
that the double-well system is unstable in linear stabil- 
ity analysis. These parameters correspond to the energy 
contour close to the homoclinic orbit and bifurcated fixed 
points (Fig. EJb)). The oscillation in the population im- 
balance has a striking resemblance to the Fig. [5] with 
the plot of the overlap in the lattice system and can be 
viewed as an analogue of the pulsating instability. 

Though we have presented different versions of the pul- 
sating instability in the two-site system [ (a) energy con- 
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tours and (b) the population imbalance ], they describe 
the same physics. The dynamically unstable system per- 
forms periodic oscillation where the system recedes far 
away from the unstable state and subsequently returns 
to this state. 

The multi site system basically shares the dynamics of 
the two-site system in a multi-dimensional phase space: 
Starting from random noise in the neighborhood of an un- 
stable steady state, the system evolves away from, and re- 
turns to, the initial state and the process repeats. These 
periodic recurrences occur in a 27V-dimensional phase 
space on the constant energy surface in full analogy with 
the two-site system. The two-site system is strictly pe- 
riodic since there is no motion out of the surface; a ID 
curve. However, in the multiwell case the dimension of 
the constant-energy surface is 2N — 1. Our pulsating in- 
stability strongly suggests that the system stays close to 
the homoclinic orbit while it evolves, but depending on 
the initial noise it still has a large state space to explore. 
But in the nonlinear multidimensional system the noise 
may cause the motion to deviate slightly. Upon looping 
around one of the stable fixed points, the multisite sys- 
tem therefore does not have to return to exactly where it 
started from. This may account for the slight variations 
in the period of the pulsations. 



V. TRUNCATED WIGNER APPROXIMATION 

In section II we have discussed the dynamics of a BEC 
within the classical mean-field theory. The GP equation 
can in general be very accurate in modeling a weakly in- 
teracting BEC. In optical lattices, however, the kinetic 
energy is represented by the hopping of atoms between 
adjacent lattice sites. This can be significantly reduced 
in deep lattices, resulting in enhanced effect of interac- 
tions and quantum fluctuations. In the following we in- 
clude quantum fluctuations in the atom dynamics using 
stochastic phase space methods. Within the truncated 
Wigner approximation (TWA) we unravel quantum dy- 
namics into individual stochastic trajectories and calcu- 
late expectation values of physical observables by ensem- 
ble averaging a large number of trajectories. 

For multi-mode dynamics TWA was introduced in non- 
linear optics in the studies of quantum fluctuations [39| . 
Details how to implement TWA in different atomic BEC 
systems may be found, e.g., in Refs. [H EH, 52] ■ in 
the TWA one neglects the third-order derivatives in the 
generalized Fokker-Planck type equation for the Wigner 
distribution function [47] . This allows us to write a non- 
linear stochastic differential equation for the Wigner dis- 
tribution ip w of the many-particle wavefunction. For a 
closed system, this equation is similar to the GP equation 
with stochastic initial conditions. 

Here we apply TWA formalism to quantum atom dy- 
namics in optical lattices. Both zero and finite tempera- 
ture nonequilirium dynamics has previously been success- 
fully studied in ID lattice systems in a number of works 



EE El 51 51 53- The effects of dynamical instabili- 
ties in lattices and TWA have been explicitly addressed 
in Refs. EE EE 53] • Since quantum fluctuations in an 
optical lattice can have a notable effect, we pay a special 
attention to evaluating the correct quantum statistical 
correlations for the initial state within the Bogoliubov 
approximation. The emphasis on quantum fluctuations 
is quite different from typical finite temperature domi- 
nated TWA approaches in higher dimensions [451 ]. 

We vary the effective ID interaction strength x/^ti f° r 
a fixed x, or a chemical potential EE- Quantum fluctu- 
ations become dominant in the limit of small atom num- 
bers and/or for strong effective ID interaction strength 
g. In the limit of Nt — > oo (for a fixed %) we recover the 
classical GP dynamics. 

For a closed system, where we ignore any dissipation 
terms, the TWA dynamics follows from the stochastic 
classical field equation, similar to GP equation, 

= - Jd&x + e.) + x k r w ■ (28) 

The difference from the GP evolution is that we generate 
a stochastic collection of the initial states and ip W 1S a 
classical Wigner representation of the full field operator. 
We evolve each stochastic realization of an initial state 
accordingly to Eq. (|28p and evaluate corresponding en- 
semble averages. Our TWA formalism is very similar to 
the one used in Refs. EE 51 53 , except that in each 
TWA realization we fix the total atom number [46| . 

A. Initial State 

In order to generate the initial state stochastically 
within TWA, we solve the quasiparticle excitation spec- 
trum using the Bogoliubov approximation. We again 
consider a stationary solution for a moving plane wave 
(f)on = e^-"^ , with u{p) = -2Jcosp + X - The lin- 
earized fluctuations around the stationary solution are 
obtained from 

i>n{t) = ^Qn&O + Sip n (t) , (29) 
such that the total number of condensate particles 

N c = (al&o), 

which is much larger than one. Analogously to our ear- 
lier classical Bogoliubov treatment, the fluctuation part, 
5ip n (t), can be written, in terms of quasiparticle opera- 
tors d q ,d q , as 

Hn(t) = ^KaV^"^*) + v;dJe-*«- a W). (30) 
i 

The operators d q ,d q obey Bose commutation relations, 
[a g ,afg^] = 6 q ,q'- Here the normal mode frequency Ct q , 
and the quasiparticle amplitudes u q , and v q are given 
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by Eqs. (fl4|) . (fT7|) . and (fl"8)) . The quasimomentum is 
denoted by 9. 

The total number of non-condensate particles in the 
Bogoliubov theory is given by 

N nc = £(K| 2 + K\ 2 ){a\a q ) + K\\ (31) 

n,q n,q 

with 

{a\a q ) = n=[e n «/ K * T (32) 

At T = 0, (a q a q ) = 0, and the non-condensate fraction 
is simply obtained from 

11. q 

In order to construct the initial state within TWA we 
replace the quantum operators (d q ,a q <') in Eq. ([30| by 
complex stochastic variables (a qi a*,) obtained by sam- 
pling the corresponding Wigner distribution function. 
Our formalism follows Ref. [42j |. except that here we fix 
the total atom number, so that in the TWA simulations 
the condensate and the non-condensate atom number 
fluctuations are related [4y|. In the Bogoliubov approx- 
imation the operators a q , a q behave as a collection of 
ideal harmonic oscillators. The Wigner function at T = 
reads [47[ 

W>„a;) = -exp[-2|cg 2 ]. (34) 

The function W(a q , a*) is a Gaussian with the width 1/2. 
Here the nonzero width mimicks the quantum noise. Due 
to the nonzero width of the vacuum modes in the Wigner 
distribution, each unoccupied phonon mode begins with 
uncorrelated Gaussian noise, distributed over the plane 
wave basis, and normalized to an average of a half particle 
per mode. This provides a seeding for scattering events in 
the dynamics, but in the end it is subtracted out from all 
normally-ordered quantum averages. For each stochastic 
realization, the number of non-condensate atoms reads 

N nc = £(KI 2 + K\ 2 )(a q a q --) + £ Kl 2 , (35) 

n.q ~ n,q 

which may fluctuate about the mean value X^nql^nl 2 - 
The ensemble average over many realizations is 
(a*a q )w = \- Since the total particle number Nt is 
conserved, the number of condensate atoms in each indi- 
vidual run is given by 

N C = N T - N nc . (36) 

Finally, we set 



in the initial state Eq. ([29]) . Note that, even though we 
consider a uniform system with a plane wave phonon ba- 
sis and uncorrelated noise in the initial phonon modes, 
the fixing of the total atom number introduces long wave- 
length correlations in the system between the condensate 
mode and the excited quasiparticle modes (46l |. 



B. Numerical Realization 

We study the non-equilibrium quantum dynamics of a 
BEC within TWA. We consider a BEC in a lattice ini- 
tially in a stable steady state which is crucial for the va- 
lidity of the TWA. At the beginning of the time evolution 
the lattice is driven to a dynamically unstable regime, 
for instance, by accelerating it through p = ir/2 or by 
modifying the atom-atom interactions. Here in our sim- 
ulations, we fix the initial velocity and change the value 
of the atom-atom interactions. We consider a small or a 
zero depletion of atoms from the condensate in the ini- 
tial state, so that the Bogoliubov approximation is valid. 
For the case of a non-interacting initial state the average 
number of non-condensate atoms is zero. We set the ini- 
tial momentum to be it. As an interacting initial state, we 
consider in all our simulations (unless otherwise stated) 
the rescaled interaction strengths, A(x/2J) = 0.284 cor- 
responding to the average non-condensate atom number 
N nc ~ 30. It should be noted that the critical value of 
A for the onset of the instability is 0.308. The number 
of lattice sites N is always taken to be 32, though larger 
lattices can also be simulated. In all the simulations we 
vary the total atom number and the interaction strength 
g oc x/Nt, so that x (and chemical potential) remains 
constant. Then the ratio x/N T represents the effective 
strength of the interactions in the system [lj| . The atom- 
atom interactions are turned up instantaneously to a de- 
sirable value so that the system evolves in the classical 
dynamically unstable regime. In all time evolutions we 
take A = 0.48. Choosing the initial state closer to the 
onset of the dynamical instability would have resulted in 
a larger depletion of atoms from the condensate, as the 
non-condensate atom number in the Bogoliubov theory 
diverges at the instability threshold (see Fig. [TO]) . 

For each individual realization of the time evolution of 
the ensemble of the Wigner distributed wave functions we 
sample the initial state according to the previous section. 
The generation of the initial state consists of replacing 
the operators (a,d^) by complex, Gaussian distributed 
variables (a, a*). We have used the Box-Mueller algo- 
rithm [34| for the sampling. As before we integrate the 
dynamical equation (Eq. ([28]) ) using the FFT split-step 
method [35| . 

C. Results 



(37) 



Since the TWA returns symmetrically ordered expec- 
tation values, instead of normally ordered ones, we need 
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FIG. 10: Number of non-condensate atoms given by Eq. f|33p 
as a function of the scaled interaction strength A within Bo- 
goliubov approximation. Note that the critical interaction 
strength for the onset of the instability is 3.08. For the inter- 
acting initial state in the simulations we consider A = 0.284 
that corresponds to Nnc ~ 30. 



to calculate the normally ordered expectation values from 
the simulation data [42], |4J| . Here we are only consider- 
ing the lowest energy band in the tight-binding approx- 
imation, so normally ordering the operator expectation 
values is straightforward. According to Ref . [42j , we have 
(here x always refers to one given site) the atom number 
in a lattice site 



n(x) = (V>*(x)V>(x)) w - ~ , 



(38) 



with the corresponding fluctuations 



An(x) = W<(r(xMx)) 2 ) w - (V>*(xMx))^ 



(39) 

The normalized phase coherence along the lattice follows 
from [H 



C(x,y) 



V n ( x ) n (y) 



(40) 



The overlap of the field amplitudes between times t = 
and t = t, which is a measure of the revival of the pulse, 
is given by 



f W (r) = (Y^r n (r)M0))w- 



(41) 



In Fig. [TT] we show a typical single-trajectory result 
for the overlaps of the state of the system with the initial 
state as a function of time for an interacting initial state 
with A = 0.284 that corresponds to the number of non- 
condensate atoms N nc ~ 30, and quasimomentum p — 
7r, for the various values of the total number of atoms 
(a) N T = 10 6 , (b) N T = 10 4 , (c) N T = 10 3 , and (d) 
Nt = 500. As stated earlier, we instantaneously turn the 



FIG. 11: Overlap of a single realization as a function of time 
sampled from the Wigner distribution for various numbers of 
particles (a) N T = 10 6 , (b) N T = 10 4 , (c) Nt = 10 3 and 
(d) Nt = 500. The initial state is interacting with the non- 
condensate atoms N nc . = 30. The parameters of the simula- 
tion are A = 0.48, N = 32 and p = n. There is no significant 
damping in the oscillation even if the non-condensate noise is 
substantial. 



interaction on to the value A = 0.48 so that the system 
evolves in the dynamically unstable regime. We vary the 
atom-atom interactions g and the atom number Nt but 
keep the value of \ (or A) fixed for each simulations. 

For a smaller total number of atoms, the non- 
condensate atom fraction in the initial state and the scat- 
tering length are larger and, consequently, quantum ef- 
fects are generally more observable. Each plot clearly 
indicates the pulsating instability without any notice- 
able damping. In each case the non-condensate parti- 
cles simply act as a vacuum noise in the system. Note 
that the pulsating period depends on the number of par- 
ticles; the smaller the number of particles the shorter the 
period of oscillation. This can be qualitatively under- 
stood in terms of the scattering events between conden- 
sate and non-condensate particles. The rate of scattering 
processes depends on the number of non-condensate par- 
ticles. Higher scattering rate leads to faster condensate 
depletion, and thus shorter period of the oscillation. 

Figure[T2]represents an ensemble average of the overlap 
f (t) sampled over 400 trajectories for the total number 
of atoms (a) N T = 10 4 , (b) N T = 10 3 , (c) N T = 500 and 
(d) Nt = 300 for interacting (full) and non-interacting 
(dash) initial state. For the non-interacting initial state 
the non-condensate atom number is zero whereas for 
the interacting initial state we again take A = 0.284 
that corresponds to the non-condensate atom number 
N nc ~ 30. For each simulations, both for interacting 
and non-interacting initial state, the time evolution is 
carried out by varying the Nt and g and instantaneously 
switching the interactions to a value such that A = 0.48. 
The other parameters of the simulations are N — 32 and 
p = tt. Though the sampling noise is still there, espe- 




2Jt 2Jt 

FIG. 12: Comparison of the ensemble average of the overlap 
of the state of the lattice sampled over 400 realizations for 
the initial number of non-condensate particles N nc = and 
N nc = 30, and for the total number of particles (a) Nt ~ 
10 4 , (b) N T = 10 3 , (c) N T = 500 and (d) N T = 300. The 
parameters of the simulations are A = 0.48, iV = 32 and 
p — 7r. There is no significant change in the nature of time 
evolution in two different initial state when the total atom 
number is large. However, the curves start deviating as the 
number of particles are reduced. 



single realization 
ensemble average 




FIG. 13: Comparison of the overlaps in a typical single re- 
alization versus an ensemble average sampled over 400 re- 
alizations with the initial interacting state for N nc = 30 
and Nt = 10 6 . The parameters of the simulations are 
A = 0.48, N = 32 and p = tt. 



cially for smaller atom numbers, the figures clearly show 
a damping in the pulsation which contrasts the time evo- 
lution in the single realization. The figure also shows that 
the smaller the atom numbers, the higher the damping 
rate. Moreover, there is no noticeable difference in the 
time dynamics for the two different initial states as long 
as the atom numbers are large. However, there is a sig- 
nificant deviation in the time dynamics for these two dif- 
ferent initial states for smaller numbers atoms. Fig. [T3] 
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FIG. 14: Ensemble average of the relative number fluctua- 
tions sampled over 400 realizations in the central lattice site 
as a function of time for the total numbers of paricles (a) 
N T = 10 4 , (b) Nt = 10 3 , (c) Nt = 500 and (d) N T = 300 
with N nc — 30. The parameters of the simulations are 
A = 0.48, TV = 32 and p = tt. 

represents a comparison of the overlaps in a typical sin- 
gle realization versus that in an ensemble average for the 
larger atom number Nt = f0 6 , and with the interacting 
initial state. The other parameters of the simulations are 
same as in Fig. (| 12|) . 

Similar sort of damping is also observed in the number 
fluctuations at a given site. In Fig. Q3]we show ensemble 
averages of the relative the number fluctuations An/n av 
at the central lattice site for the total numbers of the 
atoms (a) N T = 10 4 , (b) N T = 10 3 , (c) iV T = 500 and 
(d) Nt — 300 for an interacting initial state. All the 
parameters of the simulations including the initial state 
are the same as in Fig. (fl~2"|) . In addition to the damping 
in the oscillations the figure also reveals that the mean 
relative number fluctuations are larger as the number of 
atoms gets small. 

Our main focus in this work is to study the inher- 
ent quantum effects on the pulsation phenomenon of the 
BEC in an optical lattice. For that purpose it is obvious 
to look into the various uncertainties associated with the 
pulse. For instance we calculate the amplitude uncer- 
tainty of the first pulse in the pulsating instability as a 
function of the total number of atoms 

A(A) = yJ(A*) w (A)* w . (42) 

Fig. [15] represensts the amplitudes of the first pulse 
and the corresponding uncertainties sampled over 400 
realizations as a function of the total number of atoms, 
both for the interacting (full) and the non-interacting 
(dash) initial states. The parameters of the simulations 
are the same as in Fig. [12] For the non-interacting initial 
state the number of non-condensate atoms is zero. In 
the interacting initial state we consider the interaction 
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FIG. 15: Ensemble average of the amplitude (left) of the first 
peak and the corresponding uncertainty (right) during pulsa- 
tion as a function of the total number of atoms, Nt for in- 
teracting (full) and non-interacting (dash) initial states. The 
average is taken over 400 realizations. The parameters of the 
simulation are same as in Fig. 1121 The amplitudes remain un- 
affected for higher atom numbers, and the fluctuations acts 
as a simple random noise to initiate the instability. However, 
the quantum noise manifests itself, both in the amplitude and 
its uncertainty as the number of atoms gets small. 



parameter A = 0.284 that corresponds the number of 
non-condensate atoms N nc ~ 30. The figure shows that 
the uncertainty and the value of the pulse amplitude sat- 
urate for the higher atom numbers for both interacting 
and non-interacting initial states, but rise sharply when 
the atom number becomes small. This implies the fluc- 
tuation dominance at the smaller atom numbers in which 
the effective interaction is large. The figure also shows 
that the distinction between these two initial states is 
apparent only for smaller atom numbers. 

We see that the quantum fluctuations have an effect on 
the collapse and revival of the pulse. In the case of a sin- 
gle realization the revival seems to be very robust and re- 
peats practically forever. However, the ensemble average 
over many stochastic realizations in the Wigner method 
produces a damping in the pulsation. This damping may 
be due to the decoherence as the number of the non- 
condensate atoms grows and their interactions with the 
atoms in the condensate mode increases. The quantum 
effects are more dominant and the revival of the pulse be- 
comes progressively weaker as the number of the atoms 
becomes small. 



VI. CONCLUDING REMARKS 

In this paper we have presented new insights into the 
unstable dynamics of the BEC in an optical lattice in 
the limit of weak atom-atom interactions and by incor- 
porating quantum fluctuations. The common belief is 
that the flow of the dynamically unstable BEC in an op- 
tical lattice would be erratic, or lead to the formation 
of stable solitons. Here we moved a step further and 
show that, in the classical mean-field theory, the insta- 
bility may also trigger a quasi-periodic pulsation in the 



atom density distribution if the atom-atom interactions 
is weak. The requirement that linear stability analysis 
finds a single unstable mode gives the scale for the 'weak' 
nonlinearity and the ensuing pulsating phenomena. 

A qualitative argument has been put forward to ex- 
plain the pulsating behavior of the dynamics by com- 
paring the lattice system with the integrable double-well 
system. In the case of two wells the unstable mode leads 
to a non-trivial dynamics in the population imbalance 
such that an infinitesimal noise could produce a large- 
amplitude collective oscillation of the atoms between the 
wells. An analogous phenomenon is observed in a lattice 
in the limit of weak atom-atom interactions. We, there- 
fore, surmise that the pulsating instability is a remnant 
of the integrability. 

We incorporate the quantum fluctuations using 
stochastic phase-space methods. We use the Bogoliubov 
approximation to generate the initial state for the time 
evolution of the system. A sequence of the stochastic 
fields obtained in this way are then used to calculate the 
expectation values of the observables. We then compare 
the single realization results with the ensemble averages. 
It is observed that the quasiperiodic behavior in the time 
dynamics can still be seen in the single realizations. How- 
ever, the quantum averages show that the revival of the 
pulse becomes weaker and weaker as the atom number 
gets small. 

For experimental realizations, the flow states p ~ 7r 
near the Brillouin zone boundary can be prepared by 
accelerating the lattice [8J. Alternatively, by exploit- 
ing the symmetry of the DNLSE, every solution ip n (t) 
for the given interaction parameter A there is a solution 
(— l) n ^>*(£) for —A. That means the state for p = tt in 
the repulsive case is equivalent to the state for p = 
in the attractive case. This symmetry has already been 
used to generate solitons in a nonlattice gas (23j]. We 
speculate that the same technique can be used to ob- 
serve the pulsating instability in the lattice. However, 
given that the pulsating phenomenon only results in the 
weak nonlinearity limit, the corresponding time scales for 
the pulsation can be very long and pose a severe technical 
challenge. 

What would be seen in an experiment depends on how 
the experiment is carried out. If averaging over repeated 
experiments is called for, the ensemble averages are the 
proper quantities to compare with. On the other hand, 
it might be possible to monitor atom numbers in the 
lattice continuously, e.g., by off-resonant light scattering. 
A single realization is therefore observable as a matter of 
principles. However, the TWA scheme docs not take into 
account back-action of the measurements, which could be 
severe. This problem area will be the subject of future 
work. 
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